function fig3_12

% code to solve the BVP
%  y'' = f(x, y, y')  for  x0 < x < x1   '
%  y(x0) = y0  ,    y(x1) = y1

global ep

%ep=0.1
ep=0.01

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.05,'P',0.04)

% set boundary conditions
x0 = 0.0;	y0 = 1.0;
x1 = 1.0;	y1 = -1.0;

% parameters for calculation
nx = 600;
error = 0.000001;

% start off with a linear solution
x=linspace(x0,x1,nx+2);
gam=1;
for ix=1:nx+2
    %	y(ix)=-gam*x(ix)*(1-x(ix));
    y(ix)=1-2*x(ix);
end
dx=x(2)-x(1);
dxx = dx*dx;
err=1;

counter=0;
while err > error

    % calculate the coefficients of finite difference equation
    a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
    for j = 2:nx+1
        jj=j-1;
        z = (y(j+1) - y(j-1))/(2*dx);
        a(jj) = 2 + dxx*fy(x(j), y(j), z);
        c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z);
        v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z);
    end

    % Newton iteration
    v(1) = v(1)/a(1);
    u(1) = - (2 + c(1))/a(1);
    for j = 2:nx
        xl = a(j) - c(j)*u(j-1);
        v(j) = (v(j) - c(j)*v(j-1))/xl;
        u(j) = - (2 + c(j))/xl;
    end
    vv = v(nx);
    y(nx+1) = y(nx+1) + vv;
    err = abs(vv);
    for jj = nx:-1:2
        vv = v(jj-1) - u(jj-1)*vv;
        err = max(err, abs(vv));
        y(jj) = y(jj) + vv;
    end
    counter=counter+1;
end

newton_iterations=counter

% asy solution
for ix=1:nx+2
    dd=-0.75*(2*x(ix)-1)/ep;
    ya(ix)=x(ix)+1-3/(1+exp(dd));
end

% plot solutions
plot(x,ya,'r','LineWidth',1.2)
hold on
plot(x,y,'--b','LineWidth',1.4)
axis([0 1 -2 2])
grid on
hold off

set(gca,'ytick',[-2 -1 0 1 2])
set(gca,'YTickLabel',{'$-2$';'$-1$';'0';'1';'2'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('Solution','Interpreter','latex','FontSize',14)

say=['$\varepsilon = \,\,$',num2str(ep)];
text(0.04,-1.3,say,'FontSize',14,'FontWeight','bold','Interpreter','latex')

legend(' Composite\,',' Numerical ','Location','NorthEast','FontSize',13,'Interpreter','latex')

% if ep==0.1
%     exportgraphics(gcf,'/Users/mark/Desktop/fig29b.pdf')
% elseif ep==0.01
%     exportgraphics(gcf,'/Users/mark/Desktop/fig29a.pdf')
% end

function g=f(x,y,z)
global ep
g=(z-1)*y/ep;

function g=fy(x,y,z)
global ep
g=(z-1)/ep;

function g=fz(x,y,z)
global ep
g=y/ep;







